//
// perc2 - calculation of the fractal dimension of correlated
// bond percolation cluster hulls
//
// Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

#undef LAYERCOUNT
#define LAYERCOUNT 0
#undef CUMULATION
#define CUMULATION 1
#undef DISTPOWER
#define DISTPOWER 2
#define FORCE_HEIGHT 1.0

#include "common.cpp"

static void usage (const char *argv0)
{
  fprintf (stderr, "%s: %s [-r seed] [-R rnd.bin] [-o testcellcov.csv] [-g generator] [-s samplesize]\n", argv0, argv0);
  fprintf (stderr, "    -r seed - set the random seed\n");
  fprintf (stderr, "    -o testcellcov.csv - output data to file\n");
  fprintf (stderr, "    -g generator - use the given random generator, put 'list' to see available\n");
  fprintf (stderr, "    -s sampelcount - set the sample count\n");
  fprintf (stderr, "    -R rnd.bin - random binary bitstream data\n");
  exit(-1);
}

int main (int argc, char *argv[])
{
  const char *outfn = "testcellcov.csv";
  int rnd_seed = 0;
  const char *rndpath = 0;
  EXPONENT = 0.75;
  int samplesize = 65536;

  const gsl_rng_type * rng_type = gsl_rng_default;

  for ( int i = 1; i < argc; i++ )
    {
      if ( argv[i][0] != '-' || argv[i][1] == 0 || argv[i][2] != 0 )
          usage (argv[0]);

      switch (argv[i][1])
        {
          case 'R': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            rndpath = argv[i];
            break;
          }
          case 'o': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            outfn = argv[i];
            break;
          }
          case 'r': {
            i++;
            if ( i == argc || sscanf (argv[i], "%d", &rnd_seed) != 1 )
                usage (argv[0]);
            break;
          }
          case 's': {
            i++;
            if ( i == argc || sscanf (argv[i], "%d", &samplesize) != 1 || samplesize < 1 )
                usage (argv[0]);
            break;
          }
          case 'g': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            rng_type = 0;
            for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
              {
                if ( strcmp ((*t)->name, argv[i]) == 0 )
                  {
                    rng_type = *t;
                    break;
                  }
              }
            if ( !rng_type )
              {
                printf ("Random generators available:\n");
                for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
                    printf ("  %s\n", (*t)->name);
                exit(0);
              }
            break;
          }
          default: {
            usage (argv[0]);
            break;
          }
        }
    }

  rndf.seed (rng_type, rnd_seed, rndpath);

  printf ("Sample count: %d\n", samplesize);

  LAYERSTART = 0;
  LAYEREND = 0;

  Field *f = new Field(false);

  double start = 0;
  double end = 3;
  int steps = 60;

  std::vector<aggregator<2, true> > out;
  out.resize (steps + 1);

  for ( int i = 0; i < samplesize; i++ )
    {
      f->reset();
      double oroot = v2d(f->value(FIELD_SIZE / 2, FIELD_SIZE / 2));
      oroot += rndf.rnd() * 1;
      long double b_oroot = (oroot < 0) ? -1 : 1;
      for ( int d = 0; d <= steps; d++ )
        {
          double dist = start + d * (end - start) / steps;
          double dir = rndf.rnd1() * 2 * M_PI;
          double x = FIELD_SIZE / 2 + dist * cos(dir);
          double y = FIELD_SIZE / 2 + dist * sin(dir);
          double other = v2d(f->value(x, y));
          other += rndf.rnd() * 1;
          long double b_other = (other < 0) ? -1 : 1;
          out[d].add (oroot * other, b_oroot * b_other);
        }
    }
  FILE *fp = fopen (outfn, "w");
  for ( int d = 0; d < steps + 1; d++ )
    {
      double dist = start + d * (end - start) / steps;
      fprintf (fp, "%g,%s\n", dist, out[d].sprint());
    }
  fclose (fp);

  delete f;

  return 0;
}

